---
title: "Data Science Transport – Second Assignment – Group 12"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
library(gtfsrouter)
library(tidyverse)
library(tidytransit)
library(sf)
library(tmap)
library(units)
library(RColorBrewer)
tmap_mode("view")
```

last time {data-icon="fa-hourglass-half"}
=====================================

ADD DATA SOURCES!!!!!!!!!!!!!!!

- test1
- test2

### Chart 0
```{r}

library(leaflet)
leaflet() %>%
  addTiles() %>%
  addMarkers(lng=174.768, lat=-36.852, popup="The birthplace of R")
```

### Chart 1

```{r}
# 1. Plot the dots themselves
```

### Cars

```{r}
knitr::kable(mtcars)
```

isochrones {data-icon="fa-expand-arrows-alt"}
=====================================

```{r, include = FALSE}
##############################################################
#
#   READ GTFS DATA
#
##############################################################
# set work directions
setwd_gtfs <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/data/vbb-gtfs")}
setwd_data <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/data")}
setwd_work <- function(){setwd("~/Documents/Uni/Master/DataScienceTransport/assignment_2")}

setwd_work

# read gtfs data for monday
file <- file.path("~/Documents/Uni/Master/DataScienceTransport/data/vbb-gtfs/2020-12_2020-12-28.zip")
gtfs <- extract_gtfs(file) %>% gtfs_timetable(day = 2)

##############################################################
#
#   SET TIMES
#
##############################################################
start_time <- 7 * 3600 + 1200
end_time <- 8 * 3600

# create isochrone
# ic <- gtfs_isochrone (gtfs,
#                       from = from,
#                       start_time = start_time,
#                       end_time = end_time)

##############################################################
#
#   CREATE STOPS SF OBJECT
#
##############################################################
stops <- st_as_sf(gtfs$stops,
                   coords = c("stop_lon", "stop_lat"),
                   crs = 4326) %>% 
  st_transform(25833)

##############################################################
#
#   SHAPE DISTRICTS NEW (+ area)
#
##############################################################
setwd_data()
shape_districts_new <- read_sf(dsn = "LOR_SHP_2019-1", layer = "Planungsraum_EPSG_25833")
setwd_work()

shape_districts_new <- shape_districts_new %>% 
  group_by(BEZIRK) %>% 
  summarise() %>% 
  filter(!is.na(BEZIRK)) %>% 
  rename(NAME = BEZIRK) %>% 
  mutate(AREA = st_area(geometry)) %>% 
  select(NAME, AREA, everything()) %>% 
  mutate(AREA = (AREA / 1000000) * as_units("km2"))

# setting crs of polygons
st_crs(shape_districts_new$geometry) <- 25833

shape_berlin <- st_union(shape_districts_new)

##############################################################
#
#   SPECIFIC SHAPES AND STOPS
#
##############################################################

stops_in_berlin <- stops %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE)

# get isochrone area
# ic = gtfs_isochrone (gtfs,
#                      from = "Berlin, Sowjetisches Ehrenmal",
#                      #from_is_id = TRUE,
#                      start_time = start_time,
#                      end_time = end_time)$hull$area

##############################################################
#
#   CALCULATE ISOCHRONES
#
##############################################################
# # the following code calculates the isochrones
# # instead of running the code, we read in the pre-calculated file

# stops_ic_area <- vector(mode = "double")
# 
# # create isochrone areas for stops in 50 minutes
# for (stop_name in stops$stop_name){
# 
#   tryCatch( {
#     ic_area <- gtfs_isochrone (gtfs,
#                                from = stop_name,
#                                #from_is_id = TRUE,
#                                start_time = start_time,
#                                end_time = end_time)$hull$area
#     if(is.null(ic_area)) {
#       stops_ic_area <<- rbind(stops_ic_area, 0)
#       print(paste(stop_name, ": ", ic_area, "!!!!!!!!!!"))
#     } else {
#       stops_ic_area <<- rbind(stops_ic_area, ic_area)
#       print(paste(stop_name, ": ", ic_area))
#     }
#     },
#     error = function(e) {
#       stops_ic_area <<- rbind(stops_ic_area, 0)
#       print(paste("ERROR!!!", stop_name))
#       }
#     )
# }
# 
# ##############################################################
# #
# #   CLEANING
# #
# ##############################################################
# 
# # merge and clean
# # https://r-spatial.github.io/sf/reference/bind.html
# # https://cran.r-project.org/web/packages/units/vignettes/units.html
# rownames(stops_ic_area) <- NULL
# stops_area <- st_sf(data.frame(stops, stops_ic_area / 1000000)) %>%
#   rename(ic_area = stops_ic_area.1e.06,
#          id = stop_id,
#          name = stop_name,
#          parent = parent_station) %>% 
#   select(id, name, parent, ic_area) %>% 
#   mutate(ic_area = ic_area * as_units("km2"))
# 
# # save
# # https://r-spatial.github.io/sf/reference/st_write.html
# st_write(stops_area, "output_stops_ic_area.shp")
stops_area <- st_read("output_stops_ic_area.shp")

# more cleaning for plot
# https://dplyr.tidyverse.org/reference/distinct.html
stops_area = 
  stops_area %>% 
  select(name, ic_area) %>% 
  distinct(name, .keep_all = TRUE)

stops_area_berlin <- stops_area %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  select(-inside_berlin) %>% 
  mutate(id = paste(name, ": ", round(ic_area)))
```

Column {data-width=100}
-------------------------------------
    
### About Isochrones

WTF are you doing here?
The [gtfsrouter](https://atfutures.github.io/gtfs-router/) package allows us to calculate all stations reachable within a specified time period from a nominated station ([isochrones](https://atfutures.github.io/gtfs-router/reference/gtfs_isochrone.html)). We use the hull polygon as an indicator for the city area reachable. The time of interest is Monday from 07:20 – 08:00.


### 30min IC home trip (18:00) from Helmholtzstr.
    
```{r}
ic_einstein <- gtfs_isochrone(gtfs,
                              from = "Berlin, Helmholtzstr.",
                              start_time = 18 * 3600,
                              end_time = 18 * 3600 + 1800)

tm_basemap(leaflet::providers$OpenStreetMap.DE) +
  tm_shape(ic_einstein$hull) + 
  tm_polygons(col = "red",
              alpha = 0.2,
              border.col = "red",
              id = NA) +
  tm_shape(ic_einstein$routes) +
  tm_lines(id = NA) +
  tm_shape(ic_einstein$end_points) +
  tm_dots(col = "red", id = NA) + 
  tm_shape(ic_einstein$start_point) + 
  tm_dots(col = "green")
  
```
    
### Chart 3
    
```{r}
tm_shape(shape_districts_new) +
  tm_polygons(alpha = 0,
              popup.vars = c("area" = "AREA")) +
  tm_shape(stops_area_berlin) +
  tm_dots(col = "ic_area",
          id = "name",
          popup.vars = c("area" = "ic_area"),
          size = 0.07,
          border.lwd = 0.3,
          legend.hist = TRUE,
          n = 15,
          title = "isochrone area [km^2]") +
  tm_view(bbox = shape_berlin)
```


Column {data-width=350}
-------------------------------------
   
### area size of the hull enclosing the routed points

```{r}
##############################################################
#
#   PLOT
#
##############################################################

tm_shape(shape_districts_new) +
  tm_polygons(alpha = 0,
              popup.vars = c("area" = "AREA")) +
  tm_shape(stops_area_berlin) +
  tm_dots(col = "ic_area",
          id = "name",
          popup.vars = c("area" = "ic_area"),
          size = 0.07,
          border.lwd = 0.3,
          legend.hist = TRUE,
          n = 15,
          title = "isochrone area [km^2]") +
  tm_view(bbox = shape_berlin)
```


traveltimes to center {data-icon="fa-stopwatch"}
=====================================

```{r, include = FALSE}
##############################################################
#
#   SHAPE CENTER AREAS
#
##############################################################
# "Zentrentragender Stadtraum mit höchster / hoher Urbanität"
# of Zentrumsbereichskernen
# see page 39: https://www.stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/download/zentren/2011-07-31_StEP_Zentren3.pdf
# or page 45 (less detailed): https://www.stadtentwicklung.berlin.de/planen/stadtentwicklungsplanung/download/zentren/StEP_Zentren_2030.pdf
# recreated with QGis

shape_center <- read_sf(dsn = "shape_center_areas", layer = "center_areas") %>% 
  mutate(name = c("east", "west")) %>% 
  select(name)

shape_center_east <- shape_center %>% filter(name == "east")
shape_center_west <- shape_center %>% filter(name == "west")

##############################################################
#
#   READ GTFS DATA
#
##############################################################
# now we work with tidytransit
# calculation of shortest tt from all station to specific ones is more convinent

setwd_gtfs()
gtfs <- read_gtfs("2020-12_2020-12-28.zip")
setwd_work()

# http://tidytransit.r-transit.org/reference/filter_stop_times.html
stop_times_filtered <- filter_stop_times(gtfs, "2021-01-18", "06:00:00", "07:55:00")

##############################################################
#
#   GET STOPS
#
##############################################################
stops <- st_as_sf(gtfs$stops, coords = c("stop_lon", "stop_lat"), crs = 4326) %>%
  st_transform(25833) %>% 
  select(stop_name) %>%
  rename(name = stop_name) %>%
  distinct(name)

stops_berlin <- stops %>% 
  mutate(inside_berlin = st_within( geometry, shape_berlin )) %>% 
  mutate(inside_berlin = !is.na( as.numeric( inside_berlin ))) %>% 
  filter(inside_berlin == TRUE) %>% 
  select(name)

stops_center <- stops %>% 
  mutate(inside_center = st_within( geometry, shape_center )) %>% 
  mutate(inside_center = !is.na( as.numeric( inside_center ))) %>% 
  filter(inside_center == TRUE) %>% 
  select(name)

stops_center_east <- stops %>% 
  mutate(inside_center_east = st_within( geometry, shape_center_east )) %>% 
  mutate(inside_center_east = !is.na( as.numeric( inside_center_east ))) %>% 
  filter(inside_center_east == TRUE) %>% 
  select(name)

stops_center_west <- stops %>% 
  mutate(inside_center_west = st_within( geometry, shape_center_west )) %>% 
  mutate(inside_center_west = !is.na( as.numeric( inside_center_west ))) %>% 
  filter(inside_center_west == TRUE) %>% 
  select(name)

##############################################################
#
#   TT calculation
#
##############################################################
# what are the tt to the center areas?
# according to Nahverkehrsplan Berlin 2019-2023: ANlage 1 - Monitoringbericht (p. 12)
# standard: tt_max = 3600, n_transfer_max = 2, n_realise_stations = 0.95

tt <- travel_times(
  stop_times_filtered,
  stops_center$name,
  time_range = 5400,
  arrival = TRUE,
  max_transfers = 2,
  # max_departure_time = NULL,
  return_coords = TRUE,
  return_DT = FALSE
)

# clean it for plot
tt <- tt %>% 
  rename(from = from_stop_name,
         to = to_stop_name,
         tt = travel_time,
         departure = journey_departure_time,
         arrival = journey_arrival_time
         ) %>% 
  select(-c(from_stop_id, to_stop_id, to_stop_lat, to_stop_lon)) %>% 
  st_as_sf(coords = c("from_stop_lon", "from_stop_lat"),
           crs = 4326) %>% 
  st_transform(25833) %>% 
  mutate(tt = set_units(round(tt/60, 2), "min"))
```

Column {data-width=100}
-------------------------------------
    
### About 

WTF are you doing here?

### Chart 1
    
```{r}
```
    
### Chart 3
    
```{r}
```


Column {data-width=300}
-------------------------------------
   
### Chart 2

```{r}
##############################################################
#
#   PLOT
#
##############################################################

# https://campus.datacamp.com/courses/visualizing-geospatial-data-in-r/raster-data-and-color?ex=9
rdylgn <- rev(brewer.pal(7, "RdYlGn"))

# https://leaflet-extras.github.io/leaflet-providers/preview/
# https://tlorusso.github.io/geodata_workshop/tmap_package
# https://www.rdocumentation.org/packages/tmap/versions/3.0/topics/tm_basemap
# https://rdrr.io/cran/tmap/man/tm_view.html
# https://leafletjs.com/reference-1.3.4.html#map-methods-for-modifying-map-state

tm_basemap(leaflet::providers$CartoDB.DarkMatter) +
  tm_shape(shape_districts_new) + 
  tm_polygons(alpha = 0,
              lwd = 1.5,
              border.col = "white") +
  tm_shape(shape_center) +
  tm_polygons(alpha = 0.2,
              col = "red",
              border.col = "red"
              ) + 
  tm_shape(tt) +
  tm_dots(col = "tt",
          style = "fixed",
          breaks = c(0, 10, 20, 30, 40, 50, 60, 120),
          labels = c("0 – 10", "10 – 20", "20 – 30", "30 – 40", "40 – 50", "50 – 60", "> 60"), 
          id = "from",
          palette = rdylgn,
          title = "traveltime [min]",
          popup.vars=c("from" = "from",
                       "to" = "to", 
                       "traveltime" = "tt",
                       "departure at" = "departure",
                       "arrival at" = "arrival",
                       "number of transfers" = "transfers")) +
  tm_view(bbox = shape_center)
```


***

https://rstudio.github.io/leaflet/

- Interactive panning/zooming

- Compose maps using arbitrary combinations of map tiles, markers, polygons, lines, popups, and GeoJSON.

- Create maps right from the R console or RStudio

- Embed maps in knitr/R Markdown documents and Shiny apps

- Easily render Spatial objects from the sp package, or data frames with latitude/longitude columns

- Use map bounds and mouse events to drive Shiny logic